A real-space, real-time method for the dielectric function 
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Abstract 

We present an algorithm to calculate the linear response of periodic systems 
in the time-dependent density functional theory, using a real-space represen- 
tation of the electron wave functions and calculating the dynamics in real 
time. The real-space formulation increases the efficiency for calculating the 
interaction, and the real-time treatment decreases storage requirements and 
allows the entire frequency-dependent dielectric function to be calculated at 
once. We give as examples the dielectric functions of a simple metal, lithium, 
and an elemental insulator, diamond. 
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I. INTRODUCTION 



Real-space methods have proven their utihty in calculating the linear response of finite 
systems in time-dependent density functional theory However, there has been the 

perception that real-space methods are unsuitable for infinite periodic systems. The problem 
is that the long range polarization currents are important but are dynamically independent 
of the local state of the electrons within the unit cell. Stated differently, the polarization 
gives rise to a surface charge at the surface of any finite sample, but the resulting electric 
field is independent of the charge density within any cell in the interior. The necessity to 
introduce the polarization as an independent degree of freedom has been well recognized in 
the literature of the density functional theory 

We will show here that in fact it is straightforward to treat infinite systems in the real- 
time formulation of TDDFT, simply by adding as one additional dynamic variable the surface 
charge. Formally, this is conveniently done adding a gauge field in within a Lagrangian. A 
gauge formalism has also been very recently applied by Kootskra et al [0, however using 
a frequency representation rather than solving real-time equations. In our formulation, we 
derive the dynamic equations from the Lagrangian: 

L= I d\{ Si|V0iA-eA2:(/.,p ^ l.^v{r) ■ VV{r) + en{r)V{r) + en,Ur)V{r)+ (1) 

v.M01 + v.4P(n.')l)-^(§)^-.X<<Vi:,-f. 

Here the (pi are the Bloch wave functions of the electrons, normalized so that n{r) = 
J2i \4'i{^)\'^ is the electron density. The volume of the unit cell is Q. The electromagnetic in- 
teraction is separated into a Coulomb field V{f) that satisfies periodic boundary conditions 
in the unit cell and a vector gauge field zA(t). Note that the gauge field is uniform, without 
any dependence on r. The electric field is then given by 

S = -^v -z—. 

dt 

In these formal equations, we use units with h = c = 1. 



The other pieces of the first integral are the usual terms in the Kohn-Sham energy func- 
tional. The term en{r)V{r) gives the direct Coulomb interaction of the electrons, except for 
the surface charging. The ionic interaction is separated[| into a long-range part that can be 
associated with an ionic charge density nion{r) and a short-range part Vion- The latter de- 
pends on the orbital angular momentum of the electrons in typical ab initio pseudopotentials. 
It therefore depends on the full one-electron density matrix p{r,r') = I]i 01 We 
have emphasized this point because nonlocal interactions do not respect gauge invariance. 
The invariance is of course restored if the density matrix is gauged. If one expresses the 
nonlocality in terms of the operator V, one makes the replacement V — >■ V — ieAz. Finally, 
the Vxc is the usual exchange-correlation energy density of density functional theory. 

Requiring the Lagrangian action to be stationary gives equations of motion for 0j and 
A and the Poisson equation for V. The dynamic equation for (pi is the time-dependent 
Kohn-Sham equation, 

- 7^0i :^V^0i + —A (pi + [eV + — h ^—)(pi = iw:(Pi (2) 

im mi 2m on on at 

The equation for A is 

l^^--E(0^|V.Al0.) + -^Ar, + A / v.o„dV = (3) 
47r dt'^ m m oA Jn 

where N^. = d^rn{r) is the number of electrons per unit cell. 

II. LINEAR RESPONSE, SUM RULE AND SIMPLE MODELS 

The calculation of the dielectric function using the above real-time dynamic equations 
is very similar to the corresponding calculation of dynamic polarizability of finite systems 
0. We first solve the static equations (with A=0) to get the ground state electron orbitals 

^The separation is somewhat arbitrary, but is useful because the periodicity of V then takes the 
ionic lattice into account automatically. 
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and the periodic Coulomb potential V. The system is then perturbed by making a sudden 
change in A, A{t = 0+) = Aq. This corresponds to applying a short-duration electric field 
at t = 0, S{t) = — Ao5(t). The dynamic equations are then applied to evolve the variables 
in time, finding the time evolution of the polarization electric field S{t) = —dA{t)/dt. The 
dielectric function e{uj) is just the ratio of the Fourier components of the external and the 
total fields; it is given by 



e{uj) Aq Jo+ dt ^ ' 

Here is a small quantity to establish the imaginary part of the response. In principle the 
resulting theory automatically respects the Kramers-Kronig relation. 

The linear energy-weighted sum rule is easily derived in this formalism. The sum rule 
may be expressed as 

/ uj\m e \Ajj)dijj = — — . (5) 

Jo mil 

To calculate the sum rule with our Lagrangian, we write the integral using eq. (4) 

r culm t-\u)du = ^ r dt^lm H dcucue'^'-^' = { . (6) 
Jo ^ ' AqJo dt Jo 2Ao\dt^J^^^^ 

The second derivative in the last expression can be easily found from eq. (3). At t = 0+, 
the wave functions have not yet had time to change, ^4(0+) = Aq and (V^) = 0. Then if the 
last term in eq. (3) can be neglected, 

d^A _ A'Kc^NeAo 

and eq. (5) follows immediately. Thus the time-dependent treatment satisfies the sum rules 
automatically to the extent permitted by the last term. That term is only nonzero for nonlo- 
cal pseudopotentials, and in fact it may improve the accuracy of the theory by incorporating 
effects of the core electrons on the dynamic properties ||^. 

Let us now see how the gauge field treatment works in a simple analytically solvable 
model, namely the electron gas. As mentioned before, when the field Aq is applied, there is 



no immediate response to the operator V, since the wave function does not change instanta- 
neously. However, in the Fermi gas, the single-particle states are eigenstates of momentum 
so the response remains (V) = for all time. Putting this in eq. (3), and dropping the 
pseudopotential term, the equation for A becomes simple harmonic motion, with the solution 

A{t) = Aq cos tUpit (8) 

where Upi is the plasmon frequency. 

The dielectric function may now be calculated from the time integral eq. (4). One obtains 
the familiar electron gas result, 

eM = 1-4- (10) 

One sees that the derivation here is much simpler than the usual one using the Coulomb 
gauge. There one formulates the response in a particle-hole representation, and takes the 
external field to be of the form e*'^ '" with q finite. The dielectric function is then found by 
taking the g — limit. 

We can make another simple model for the opposite extreme of a tightly bound electron 
in the unit cell. Assume that the ion potential eVion{r) + 5Vion/5n can be approximated 
by a harmonic oscillator potential in the region over which the electron wave function is 
appreciable. According to Kohn's theorem the response is just the same as for an 
isolated electron in the same ionic potential. This comes out of eq. (2-3) in the following 
way. The initial impulse Aq starts the electron moving, and as a result both V{r) and 
^Vi^c/^^l'") become time-dependent. Together with the changing A, the electron in the unit 
cell drags its self-induced field with it, and the accelerations associated with these three 
terms in eq. (2) exactly cancel. The remaining ionic terms then produce simple harmonic 
motion for ch. 
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III. NUMERICAL DETAILS 



The computational algorithm we employ is identical to the ones we used for clusters and 



molecules, which was based on a method introduced in nuclear physics |10]. The Kohn 



Sham operator is represented on a real space grid as in ref. |rT||. There are a number of 



technical details associated with the periodicity and with the gauge potential that did not 
arise for the finite-system calculations. In the new code, the potential V{r) is calculated by 
Fourier transformation of the Poisson equation rather than a relaxation method. This gives 
automatically the required periodicity to V{r). The wave functions 0j represent Bloch states 
of the periodic lattice, and they are constructed with the corresponding periodic boundary 
conditions labeled by the Bloch momenta k. The periodic boundary condition on the Bloch 
wave function ^^(r + a) = exp(iaA;)0fc(r) is easily implemented in the relaxation method used 
to find eigenstates. In practice, many Bloch states are needed to obtain smooth dielectric 
functions. However, constructing the states takes much less time than for the same number 
of electrons in a finite system, because the Bloch states in a given band are automatically 
orthogonal. 

We use here the same energy density functional that we used previously for finite sys- 
tems. Only the valence electrons are included explicitly; core electrons are treated by a 
pseudopotential [|T2],|T^. The exchange-correlation energy of the electrons is calculated in 



the local density approximation following the prescription of ref. |IJ . 

The presence of a vector gauge potential requires a modification in the pseudopotential 
calculation, as indicated in the introduction. In particular, the A-dependence of the Vion 
term in eq. (2) must be consistent with the last term in eq. (3) in order to have the 
algorithm conserve energy. We implement the A-dependence of Vion simply by gauging the 
density matrix directly, 

V^onApir,r)) = V.on(e^^(^-^'V(r,r')) 
As in the finite systems calculations, energy is conserved to a very high accuracy with the 
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algorithm [|T0| , provided the time step is less than the inverse energy span of the Kohn-Sham 
operator. 



IV. LITHIUM 

In this section we demonstrate the feasibility of the method with lithium as an example 
of a simple metal. As other alkali metals, lithium has a Fermi surface which is nearly 
spherical. However, unlike sodium and potassium, the effective mass of the electrons at the 
Fermi surface is significantly enhanced over the free- Fermi gas value (m* ^ l.Qnie)- 

The Kohn-Sham operator is represented in coordinate space with a uniform spatial mesh. 
The lattice spacing of the bcc unit cell of Li metal is a = 3.49 A, and we take a mesh spacing 
of Ax = 0.58 to subdivide the cell into a 6^ lattice of mesh points. We use a time step of 
At = 0.01 eV~^ which is sufficient to conserve energy to 10^^ eV over the time integration 
interval, T = 18 eV""*^. 

We sample the occupied states with a uniform mesh in momentum space. With a finite 
set of occupied orbitals, the allowed excitation energies will be discrete, and the metallic 
behavior, e{uj) oo, is only reached in the limit of an infinitely dense momentum space 
lattice. However, it is our view that the TDLDA loses validity at long times when other 
degrees of freedom can be excited. This is the case for the low frequency response of metals, 
where the imaginary part of the dielectric responses is dominated by phonons and inelastic 
electron scattering. 

In Fig. 1 we show the normalized time-dependent polarization field, AQ^dA/dt over the 
time interval t = [0, 18] eV~^. The inset shows an expanded view of the initial response in 
the interval [0,1] eV~^. The dashed line is the comparison with the linear behavior deduced 
from the sum rule, A^^dA/dt = —Aire^ Net/mQ. The agreement shows that the local sum 
rule in nearly satisfied, despite the fairly large optical effective mass. 

In Fig. 2 we show the inverse dielectric function computed from eq. (4) for various 
meshes in the Brillouin zone. We employ r] = 0.2 eV to smooth the response in the Fourier 
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transformation. We see that the response becomes smoother, the more finely the Fermi sea 
is sampled. With a 32^ lattice of Bloch states, we get results smooth enough to be compared 
with measurement. 

In Fig. 3 we show the real and imaginary parts of the inverse dielectric function in the 
frequency interval — 20 eV/h. The dashed lines show the empirical function from ref. W5 



There is also another theoretical calculation in the literature, |T6|. The agreement is quite 
good, especially considering that the calculation is ab initio with a energy density functional 
that much simpler than more recent ones. The main feature in the absorptive response is 
the plasmon at 7 eV and its width. The peak position is significantly downshift from the 
naive plasmon frequency, u = ^Ane'^n/m ^ 8 eV. The width is associated with interband 
transitions and is also well reproduced. 



V. DIAMOND 

In this section we compute the dielectric response of a typical elemental insulator, di- 
amond. The diamond lattice is represented in our calculation by the primitive unit cell 
which contains 8 carbon atoms. The four valence electrons of each carbon are calculated 
explicitly while the core electrons are only treated implicitly by the pseudopotential. We 
found in earlier studies of carbon structures that the Kohn-Sham Hamiltonian requires a 
mesh spacing Ax = 0.3 A to get orbital energies to an accuracy of 0.1 eV; in the calculation 
here we take a 12^ lattice in the unit cell which implies Ax = 3.56/12 A = 0.297 A. With a 
smaller mesh spacing than for lithium, the span the Kohn-Sham operator is increased and 
the time step At must be reduced accordingly. We use here At = 0.002 eV~^. 

With a cubical unit cell and 8 carbon atoms, there are 8 ■ 4/2 = 16 occupied bandsQ. In 
each band we take a lattice of up to 16^^ points to represent the Bloch states. 



^The bands are actually two-fold degenerate because we have not exploited the symmetry that 
allows a smaller unit cell with 4 carbons. 
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For an insulator, a reference point of the vector potential A{t) should be irrelevant. Since 
all the points within Brillouin zone are occupied, there is a static solution of eq.(2-3) for 
each constant vector potential. However, in our real-space, real-time calculation, we found 
this is violated due to the discrete representation in coordinate and momentum spaces. The 
total energy does depend on A, and there appears a spurious low frequency mode in the 
vector potential A{t) associated with the adiabatic evolution of the electronic wave function. 
Below that frequency, a spurious conduction feature appears in the response. 

These features are shown in the plots of the response in Fig. 4. We see that the spurious 
adiabatic evolution gives an unphysical plasmon at ^ 1.2 eV, which dominates the dielectric 
response at lower frequencies. Though the amount of strength associated with this spurious 
plasmon is very small, 0.007 electrons out of the total of 32, it has a qualitative effect on the 
response at very low frequency. The frequency and strength decrease the finer the spatial 
mesh, showing that it is an artifact of the discrete mesh representation of the coordinate 
space. 

To infer the dielectric function near u = 0, we apply the Kramers-Kronig relation to the 
imaginary part of the response, but excluding the spurious plasmon peak. This gives the 
predicted dielectric function shown in Fig. 5. The empirical dielectric function is shown as 
the dashed line. The agreement is good, as indeed was found solving the TDLDA equations 
by other methods [jl^, but one can also see the effect of the well-known shortcoming of 



TDLDA, that the predicted band gaps are too small [|T8|,[T9|. The theoretical absorption 
strength become significant starting at about 5 eV excitation, while the empirical absorption 
begins at around 7 eV. Nevertheless, the dielectric constant comes out in good agreement 



with the empirical [|T7|, being within a percent of the empirical value of e(0) = 5.67. 



VI. CONCLUSIONS 

We see that the method not only works in principle, but produces fairly accurate di- 
electric functions in the cases of a simple metal and a simple insulator. In lithium, the 
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theory describes the metaUicity as well as the interband transitions. In diamond, there is 
a spurious plasmon at low frequency due to the discrete mesh representation in coordinate 
and momentum spaces. However, it can be easily dealt with and then the dielectric func- 
tion has an excellent quality except for a small band gap region. We find that two benefits 
of the real-space, real-time formulation of the TDLDA in finite systems [^] are preserved 
in our implementation here. The real-space method allows the Kohn-Sham operator and 
electron-electron interactions to be evaluated efficiently [|ll|. Computational efficiency is 
also gained by calculating the response in real time in that all frequencies are calculated 
at once. Finally, the method requires much less storage than methods using a particle-hole 
representation of the time- varying wave function. 
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FIGURES 
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FIG. 1. The induced polarization field {l/AQ)dA/dt in lithium metal is shown as a function 
of time (solid). In this calculation, the occupied states were represented by points on a 16^ mesh 
in A;-space. The dashed line show the early-time behavior required by the sum rule. 
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FIG. 2. Real and imaginary parts of the inverse dielectric function in lithium, shown for 
various meshes on the Brillouin zone. 
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FIG. 3. Real and imaginary parts of the response e ^{co) as a function of frequency. Here 
the orbitals of the valence band were presented by Bloch states on a 32^ mesh in A;-space. The 



empirical response from ref. |15] is shown with the dashed lines. 
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FIG. 4. Real and imaginary parts of the response e ^{u) for diamond. Here the orbitals of 

the valence band were represented by Bloch states on a 16^ mesh in /c-space. 
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FIG. 5. Real and imaginary parts of the dielectric function e{Lo) for diamond. The spurious 
plasmon has been excluded by using the Kramers-Kronig relation to determine the real part of the 
dielectric function, integrating over the imaginary response from 4 eV. The dashed curve shows 
the experimental dielectric function, taken from ref. |15|. 
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